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Abstract 

In turbulent flows, local velocity differences often obey a cascade-like hierarchical 
dynamics, in the sense that local velocity differences at a given scale k are driven by 
deterministic and random forces from the next- higher scale k — 1. Here we consider 
such a hierarchically coupled model with periodic boundary conditions, and show that 
it leads to an iV-th order initial value problem, where N is the number of cascade steps. 
We deal in detail with the case N = 7 and introduce a non-polynomial spline method 
that solves the problem for arbitrary driving forces. Several examples of driving forces 
are considered, and estimates of the numerical precision of our method are given. We 
show how to optimize the numerical method to obtain a truncation error of order 0(h 5 ) 
rather than 0(h 2 ), where h is the discretization step. 

1 Introduction 

Hierarchical dynamics arise quite commonly for complex systems that consist of many 
subdynamics that drive each other in a selfsimilar way. Often, the dynamics at a given 
scale k couples to the dynamics at a higher scale k — 1 in a simple way, and these types of 
problems can sometimes be dealt with in an analytic way. A typical example are cascade 
models in turbulence (see, e.g. [USE] and references therein), which are a useful tool to 
characterize the selfsimilar features of turbulent flows at high Reynolds numbers. 

In this paper we consider a very simple hierarchical model that can be physically inter- 
preted as representing velocity differences at different scales that are driven by deterministic 
and stochastic forces from the next higher scale, in a medium with viscosity and a given 
driving force at the top scale. We show that the problem reduces to the solution of an initial 
value problem associated with a differential equation of order N, where iV is the number 
of cascade steps. As an example, we deal in detail with the case N = 7, which corresponds 
to typical values of cascade sizes observed in turbulent flows. There are three standard ap- 
proaches to solve initial value problems numerically, the finite difference method, the finite 
element method and the spline approximation methods. We introduce a non-polynomial 
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spline method for the numerical solution of this initial value problem, and estimate the 
precision of our numerical treatment. Higher-order boundary value problems are effectively 
solved using non-polynomial spline methods El El E] . It turns out there is an optimum 
choice of the interpolation parameters where our methods yields best possible results (5th 
order rather than 2nd order in the discretization step). Some analytically solvable examples 
of driving forces are dealt with as examples, for illustration of our general method. 



2 Cascade model 

Let us consider a model of damped particles in a viscous medium that are driven by rapidly 
fluctuating forces. Suppose there are N such particles. We denote the velocity of each 
particle as y^ k \t), k = 1, 2, . . . , N. A very simple, uncoupled model would be that the 
particles are damped by a linear friction force, which is proportional to velocity, and a 
rapidly fluctuating driving force (t) that is independent of velocity: 



y 



(*) _ _r„( fc ) 



Ty^+L^(t) (2.1) 

r > denotes the friction constant. If L^ k \t) is Gaussian white noise, then this model just 
leads to the Ornstein-Uhlenbeck process, performed independently by each particle [4j. 

Here, however, we want to modify this model to a more interesting interacting dynamics. 
First of all, we allow L^ k \t) to be any time-dependent driving force, and assume that it 
is differentiable. Next, we construct a coupled hierarchical model, by replacing y^ on the 
right-hand side of the above equation by the nearest neighbour We may physically 

interprete y^ as a local velocity difference in a turbulent flow at spatial scale r = 2~ fc . 
The physical interpretation is that the change of velocity due to friction forces at a given 
scale k is proportional to the velocity at the next-smaller scale k + 1. This reminds us of 
the fact that in cascade models of turbulence energy dissipates from larger scales down to 
smaller scales. In our hierarchical model the actual dissipation force at spatial level k is 
proportional to the velocity difference at the next smaller scale k + 1: 

y( fc ) = -iV + (t) k = 1, . . . , N (2.2) 

We now show that this model, with a cascade of size N, leads to an initial value problem 
of the form 

^-y=(-l)»T N y + g(t), (2.3) 

where y(t) = y^'(t) is the velocity difference at the top of the cascade (k = 1), g(t) is a 
driving force at the top of the cascade, and we have assumed periodic boundary conditions 
at the top and bottom of the cascade, i.e. y^ n+l \t) = y^ l \t). 
To derive eq. (12. 3j) . we differentiate eq. (12.21) to obtain 

y< k ) = -Ty { - k+l ) +L^ k \t) (2.4) 

= +rV fc+2) -TL^ k+1 \t) +U k \t), (2.5) 



where in the last step we used eq. (|2.2p with k replaced by k + 1. Further differentiation 
yields 

&_ y (k) = r y fc+2 )+L( fc )(i)-rL( fc+1 )(i) (2.6) 
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= -ry fc+3 > + T 2 L^ k+2 \t) - TL^ k+1 \t) + L {k \t). (2.7) 



Finally, for a cascade with N steps one arrives at 

aN 

y<y {k) {t) = (-irr N y^ + g(t), (2.8) 
where g(t) is a linear combination of derivatives of the driving forces at the various scales: 

g(t) = (-rf-W+t-Vit) + (_ r) JV-2^(fc+7V- 2 ) + _ _ _ (29) 

Our implemented periodic boundary condition of the cascade y^ (t) = (t) simply means 
that at smallest scales the dynamics should just be the same as at the largest scales, which 
is a self-similarity assumption. For N odd and defining y^\t) = y(t) we arrive at 

aN 

' T y + T N y=g{t), (2.10) 



dt N! 

which has to be supplemented by a set of N initial conditions, corresponding to initial 
velocities at the various scales. This initial value problem will be solved in the next section. 



3 The initial value problem 

In the following, we choose as an example N = 7. In turbulence simulations, the driving 
forces at various scales of the turbulent flow are only known in numerical form, and given 
the chaotic nature of the forces it is important to implement high precision numerical meth- 
ods that optimize the numerical solution of driven velocity fields within a given cascading 
subdynamics. In the following, we introduce a high-precision non-polynomial spline method 
for the solution of the initial value problem (IVP) that corresponds to our cascade model. 



Consider the following seventh order initial value problem 



y (7) (t) + f(t)y(t) 



vw 
y {2 \a) 

y {A \a) 

.( 6 ) (a) 



git), 



y( 3 ) (a) 
y (5) («) 



G [a, 
= ui, 

U 5 , 



b], 



(3.1) 



where uf, i = 0, 1, . . . , 6 are finite real constants while the functions f(t) and git) are con- 
tinuous on [a,b]. This is just the dynamics derived in the previous section, provided the 
effective friction constant F depends on time t. The notation is slightly different since now 
denotes the i-ih derivative with respect to t. For our physical application in terms of a 
cascade- like model, we need to keep the functions f(t) and g{t) quite general since they are 
unknown in a turbulent flow. Our aim is to provide a proper numerical method to provide 
a most accurate solution of this IVP. It turns out that a non-polynomial spline method 
is very useful in this context. For particular choices of the interpolation parameters our 
method provides optimum results (error term of order five rather than 2), as shown in the 
following sections. 
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4 Nonpolynomial Spline Method 

To develop the spline approximation to the problem (|3.ip . the interval [a, b] is divided into n 
equal subintervals, using the grid points £j = a + ih ; i = 0, 1, . . . , n, where h = (b — a)/n. 
Consider the following restriction Si of the approximate solution S to each subinterval 
[U, t i+ x), i = 0,1,... ,n- 1, 



= aj cosw(t - tj) + 5j sinw(t - tj) + Cj(t - tj) 6 + - tj) 5 + ej(t - tj 



(4.1) 



Let 



2/i = Si(ti) 
Mi = sf (i*), 
Ui = 



mi 

N; 



0,1,..., n. 



(4.2) 



Following [6J and postulating continuous derivatives at knots, consistency relations between 
the values of splines and their seventh order derivatives at knots are obtained as 



{ah 7 U^ 7 + ph 7 Ui- & + jh 7 Ui_ 5 + 5h 7 Ui. 4 + 5h 7 Ui- 3 + -fh 7 U^ 2 + fih 7 XJ^ x + ah 7 Ui) 
[-120^„ 7 + 840 yi _ 6 - 2520 yi -5 + 4200^-4 - 4200y;_ 3 + 2520y;_ 2 - 840y,_i + 120^] ; 



7,8, 



(4.3) 



where 



a 



+ 



P 



/ 120(cos6>- 1 i (ill 

V e 7 sm6 
/ 600(1 - cos I 

V 7 sin0 



+ 



1 



5 sin 3 sin 9 69 sin 
60(2 cos 8- 3) 5(2 cos 9 



5 sin 9 



+ 



9 s 



sm I 



(2 cos 8 - 57) 
66> sin^ 



/1O8O(cos0 - 1) 180(2 cos 9 + 1) 45(2 cos + 1) (38 cos - 101) 



7 sin0 



+ 



5 sin0 



3 sin 9 



29 sin 9 



and 



/ 600(1 -cos 6) 6O(4cos0 + l) 5(20 cos 8 - 1) (604 cos 9 - 359) \ 
V 8 7 sin 8 5 sin 8 3 sin 8 68 sin 8 ) 



Here 8 = uoh is an arbitrary parameter. The relation (|4.3p forms a system of (n — 6) linear 
equations in the in) unknowns (jji, i = 1, 2, n), while Ui is taken from IVP (|3.ip to be 
equal to -fiVi + gi, i = 0, 1, . . . , n. 

Six further equations (end conditions) are required to obtain a complete solution for the yiS 
appearing in eq. (|4.3|) . These equations are calculated using the method of undetermined 
coefficients [5] as follows 



U - WUi + U 4 



30666 
8867 



u 2 + u 5 



1 

h 7 



1 

h 7 



512540 



27 



-y - 20160yi + 1260?/2 



2240 161000, (i) 
H 7, — h V 



27 



9 



-957600 



(4.4) 



8867 



1048320 
Vl + ^867- y2 



90720 
8867 



1/3 



866880 (i) 
^86T Vo 



4 



u 2 



278026 . 
94221 



U 3 + U 6 = 



U 4 + U s 



1 

it 



1209600 2 ( 2 ) 

87864 ^5(5) 
8867 

-67340 



^867~ ftyo 



352800 4 (4 ) 
^867~ ^ 2/0 



(4.5) 



80640 

2^2 + 7777777^3 



^3 + ^7 = ^ 



10469 10469" 
_4200 2 (2) 

361 Vo 
95536 ^ 5 ( 5) 

31407 Vo 

10808537040 

487056529 

8243655840 ( i) 
''//o 



700 
551 



2/4 



54040 
10469 



(i) 



20720 ,3(3) 

1653 y ° 



85400 4 (4) 
10469 Vo 



(4.6) 



1 



-2/3 



13373418240 2564881200 
-2/4 + 



487056529 

26287914240 , 2 (2 ) 
/' l/o 



487056529 



2/5 



487056529 " u 487056529 
40576352880^3 (3) 39377200800 ^ 4 (4 ) 

+ 487056529 V ° + 487056529 V ° 
25438766892 - (5) 9474762304 , „ (6 ) 



487056529 
1 I" 2645350155 
W 436783036 



487056529 
869117760 



(4.7) 



■2/4 



109195759 



2/5 + 



831120885 
436783036 



2/6 



31279815, (i) 3666455415 l2 f2 ) 
H fy/n H ^ Vn 



7530742 

3572264955 , 
+ h'yl 



218391518 
3 ( 3) 8717751945 



109195759 
3525702999 5 (5) 

109195759 



218391518 
1634628387 ^ 6 ^( 6 ) 

109195759 



(4i 



and 



U 5 + U 9 



1 



132838307280 
61865369749 



■2/5 



183300929280 
61865369749 



■2/6 



50462622000 
61865369749 



2/7 



82375685280, m 402603647040 , 2 (2 ) 
H fy/n H h 2 yX> 



61865369749 
946588828080 



61865369749 



^2/^ + 



1390554453120 



4„,(4) 



61865369749 
1350858565644 ^ , 

61865369749 



i4 V 

61865369749 Vo 
749461929944 ^ (6) 

61865369749 



(4.9) 



Basically, one does a power expansion in h, and postulates that the low orders in h vanish. 
Our calculations were done using Mathematica. 

The local truncation errors associated with the linear equations (|4.4p — (|4.9p and (|4.3p 
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are calculated as 



-5.778/iy 9 )(ti) + 0(/i 10 ), 
-6A72h 9 y^(t 2 ) + 0(h w ), 
-7.2Wh 9 y( g \t 3 ) + 0(h w ), 
-19.288/i 9 y( 9 )(t 4 ) + 0(h w ), 
-25.620k 9 yW(t 5 ) + 0(h w ), 
-33.020/iV 9 )(t 6 ) + 0(h w ), 
§(-100 + 25a + 13/3 + 5 7 + 5)h 9 y^\ti) + 0(h 



10 \ 



1, 
2. 
3, 
4, 
5, 
6, 
7, 



(4.10) 



,n 



and 



|f || = c/i 9 i? = 0(/i 9 



i? 



max \y 

te[a, b] 



(9) 



(4.11) 



where c is a constant which depends only on the values of a, f3, 7 and 5 and is independent 
of h. Moreover, a, /3, 7 and 5 are taken such that a + /3 + 7 + <5 = 60. 
In general, the solution of the system of linear equations (I4.4f) — (j4.9j) and (I4.3P is second 
order convergent. 

The local truncation error of the system (I4.3P can be expressed in the following form 



2(-60 + a + p + 7 + 6)h 7 y( 7 \ti) + (-60 + a + /3 + 7 + 5)h 8 yW(ti) 
+i(-100 + 25a + 13/3 + 57 + 5)h?y^ (U) 
+£(-120 + 37a + 19/3 + 7 7 + f)h w yW(U) 
+^(-228 + 337a + 97/3 + 17 7 + 5)h u y^(ti) 

+m(- 380 + 781a + + + <^ 12 y (12) (^) 

+0(/i 13 ), 

i = 7,8, ... ,n, 



(4.12) 



therefore, the order of the truncation error ij can be improved to be of order h , if a = 
— |, /3 = — + (5, 7 = ijSi — ^. Correspondingly, the end conditions with local 



truncation error of 0(h ) can be determined as 



1 

h7 



Ux + 
1 



24407 TT 59362 10662 TT 9 7 rT 

£/ 2 ^3 U A + f7 5 

- 109 109 



109 109 
66633336 19958400 9979200 

-yo T7^ — yi h — tt^ — 2/2 



+- 



545 

9114336 



109 



109 



4435200 1247400 

"2/3 H 77^ — 2/4 



109 



109 

202055040421 



1995840 
109 



— h\ {7) 
109 



77878525838 18661788874 
U2 + 777777=7777777 ^3 + 777777=7777777 ^4 



109 



2434662535 . 



798336 
545 2/5 

(4.13) 



554613069 184871023 184871023 554613069 

-4011644165760 8861887188000 73815832428800 

"2/1 H -.n.n^^o — 2/2 7777^777777777^; — 2/3 



■C/ 5 + C/ 6 



184871023 

4618760731200 

H "V4 — 

184871023 

4345911046400 (1) 
554613069 2/0 



184871023 1663839207 

1470208924800 1826678971040 

2/5 H 1ffMimm 2/6 



184871023 
1035868310400 

184871023 



1663839207 
:2 (2 ) 183470425600 

184871023 



h z y K 



(4.14) 



6 



13173366154319505819 TT 5923535526089565973 TT 

n^An^nnn^An^A ^3 o nn A n o A or^o n ^4 



604803696004634 
2639790737228529743 



U 5 + U 7 



1 

h7 



604803696004634 
25352931909798309915 



U 3 
1 



i 



302401848002317 



198185856313975120000 



43200264000331 2/2 " 
72535593878062755750 



2/3 



129600792000993 
44672825515677652800 



43200264000331 
44958164899589796925 



2/4 



43200264000331 
2139803134054971840 



-1/5 



-ye 



129600792000993 " 43200264000331 
5764036699720950200 , m 7374675959642702700 



2/7 



43200264000331 
3273172503578299200 , , (3) 



43200264000331 
3169885805313999875741 

17618985121607404312 
5366584014500607349360 



1,2 (2) 

43200264000331 y ° 
521467194925746900 , 4 u) 
— — h //,', 



43200264000331 
544760103861334609083 



(4.15) 



U 6 + U 7 



19821358261808329851 ~ V3 ' 
5581229831865256388400 



17618985121607404312 
1627194491533397967735555 



2/4 



+ 



+ 



2202373140200925539 
1793040053953186573920 



2/5 



ye 



2202373140200925539 
78317811652764229087465 



yi 



1127615047782873875968 
_ 324347788172453021646845 

~ 158570866094466638808 
147257531382013448691645 

~ 1127615047782873875968 
(1) 42412130734120984954425 



"2/8 



2 (2) 

845711285837155406976 " u 140951880972859234496 Vo 
13509365844145615609375,3 (3 ) 7826257773953806554675 u4 (4 ) 



35237970243214808624 
447768854035545682017 5 (5 ) 



35237970243214808624 



8809492560803702156 
6169811365491003355386625 



u 7 + u 9 



(4.16) 



364845537886699795641421 
413182203198678792199193360481 



-2/4 



1011762526223941981900336800 



-2/5 



373601830795980590736815104 " 364845537886699795641421 
998387082478934194463004566965 _ 566429213407879867786917120 

354629862825872201363461212 m 364845537886699795641421 Vl 
173230355267937275186019127455 5254626822196644195075230752 

2/8 ovnr^Anr^nAnonmnAnonr^n 2/9 



373601830795980590736815104 
1828802733123508354716945025585 ( i) 

7565437073618606962420505856 Vo 
54937023892836663800655898465 2 (2 ) 

74170951702143205513926528 V ° 
319183920456421230207708911935 



88657465706468050340865303 



-h 3 v {3) 

315226544734108623434187744 y ° 
81879551659139637198985554365 ^ (4) 

1050755 1491 1369541 144729248 Vo 
2991554077139003376141526763 5 (5) 

8756292909280795095394104 V ° 



7 



303485670688565607390252013 ^ 6 (6) " 
+ 4378146454640397547697052 Vo 



(4.17) 



and 



U 5 + U w 



1 

W 



19038680213948167651954555270266 3612553213748861716357961962885 

"2/5 n^ooo^o^oo-,^-, r^-,^^-,^ 2/6 



43087137994818537402205515625 * 2680364680381671574716400716 
4872382360659412438888663650 3637650046079073899112436800 

+ 2757576831668386393741153 Vl 2757576831668386393741153 V8 
396380487483749310137552942650 

+ 670091170095417893679100179 V9 
25866940054548307411938055846689 

172348551979274149608822062500 Vl ° 
38692597856846972600037888063421 (1) 

+ 698011635516060305915729353125 Vo 
9138473244848295227157642242648 2 (2 ) 

46534109034404020394381956875 V ° 
2879180812841847947352594547754 ^ (3) 

+ 9306821806880804078876391375 V ° 



4.,(4) 



1170428496335992779989052992 } 

+ 4278998531899220266150065 Vo 

2186365487813281315497916274 5 (5) 

+ 15909097105779152271583575 V ° 

1319003601532979667100927096 ^ 6 (6 ) 

+ 41363652475025795906117295 V ° 



(4.18) 



It turns out that if a, (3, 7 and 5 are chosen as a = W — |, /3 = — ^jjl + 5, 7 = — 
then the order of truncation error of eq. (|4.3p is 0(h 12 ) and the order of convergence can 
be improved up to five based on the improved order of the end conditions. 

To illustrate the implementation of the method, three examples are discussed in the fol- 
lowing section. We choose a driving force at the top of the cascade for which the IVP can 
be solved analytically, and then investigate the error terms of our numerical method by 
comparing with the exact solution. 
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5 Numerical Examples 



Example 1 

Consider the following initial value problem 

V (7) (t) + y(t) = -{t 2 -43) cos(t) + (-1 + t 2 -Ut)sin(t), 

y(-i) =o, 

yW(-l) = 2sm(l), 

y( 2 )(-l) = -4cos(l) - 2sm(l) 

?/ 3 )(— 1) = 6cos(l) — 6sm(l) 

y ( 4 )(-l) = 8co«(l) + 12am(l) 

y( 5 )(-l) = -20co«(l) + 10sm(l) 

y (6) (-l) = -12cos(l) - 30a*n(l). 

This basically corresponds to a periodic forcing at the top of the cascade whose amplitude 
is a particular quadratic function of t. 

The analytic solution of the above problem is 

y{t) = (t 2 - 1) sin(t) , 

meaning the velocity of the driven particle oscillates with increasing amplitude. Since we 
can solve this example analytically, we can easily determine the numerical error of our 
method. The observed maximum errors (in absolute values) associated with yi, for the 
problem (|5.1|) . corresponding to the different values of a, (3, 7 and 5, are summarized in 
Table 1. 

The observed maximum errors (in absolute values) associated with yi , for the problem ()5.1|) , 



Table 1: Maximum absolute errors for problem (|5.ip in j/j. 



n 


a 2 , /3 2 

49 £ = 51 
1 2 ' 2 


a = 0, /3 = 
7 = 0, <5 = 60 


a = 10, /3 = 10 
7 = 10, 5 = 30 


12 


2.88 x 


3.04 x 10" 1 


2.76 x 10" 1 


24 


3.09 x 10" 2 


3.56 x 10" 2 


2.73 x 10~ 2 


48 


2.5 x 10~ 3 


3.9 x 10" 3 


1.4 x 10~ 3 


96 


1.70 x 10" 4 


7.37 x 10~ 4 


3.19 x 10" 4 



corresponding to the improved end conditions are summarized in Table 2. A significant 
improvement is obtained. 



te [-1, 1], 1 



(5.1) 
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Table 2: Maximum absolute errors for problem (|5.1|) in y^. 



n 


\y(U) -yi\ 


10 


2.25 x 10" 1 


20 


2.08 x 10" 6 


40 


7.50 x 10~ 7 



Example 2 



Consider now an exponential forcing at the top of the cascade, namely the following IVP: 



(5.2) 



y (7) (t) 


-y(t) = 


-7e*(5 + 2i), 0<t<l 


1/(0) 


= o, 


y«(0) =1, 


y (2) (o) 


= o, 


2/ (3) (0) =-3 


y (4) (o) 


= -8, 


y(5)( ) =-15, 


y( 6 )(0) 


= -24. 





The corresponding analytic solution is 

y(t) 



t(l - t) e* 



In this case, for a limited amount of time, velocity differences grow exponentially, consistent 
with the chaotic nature of the velocity field in turbulent flows. In this case the observed 
maximum errors of our method are shown in table 3. 

The observed maximum errors (in absolute values) associated with m, for the problem 
(|5.2p . corresponding to the improved end conditions are summarized in Table 4. 



Example 3 

Consider the following problem: 



2/(0) 

y (2) (0) 
y {4) (0) 

yW(0) 



-- 7(-6e* + e*(l - t)) + (-7e* + e'(l - t))t, 
0, y( 1 )(0) =1, 

0, y( 3 )(0) =-3 

-8, y( 5 )(0) =-15, 

-24. 



< t < 1 



(5.3) 



The corresponding analytic solution is again 

y(t) = t{\ 



t) e* 
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Table 3: Maximum absolute errors for problem (|5.2p in y^. 



n 


a 2 , /3 2 

49 § _ 51 
/ 2 ' 2 


a = 0, /3 = 
7 = 0, 5 = 60 


a = 10, /3 = 10 
7 = 10, 5 = 30 


10 


1.5 x 10~ 3 


1.6 x 10~ 3 


1.5 x 10~ 3 


20 


1.75 x 10" 4 


1.94 x 10" 4 


1.60 x 10" 4 


40 


1.81 x 10" 5 


2.62 x 10" 5 


1.32 x 10~ 5 



Table 4: Maximum absolute errors for problem (|5.2p in y^. 



n 


\y(U) -yi\ 


10 


1.82 x 10" 1 


12 


2.15 x 10" 8 


15 


3.65 x 10" 9 



Table 5: Maximum absolute errors for problem (15. 3h in yi. 



n 


a — 2 , P — 2 

49 g _ 51 
7 2 ' 2 


a = 0, /3 = 
7 = 0, 5 = 60 


a = 10, /3 = 10 
7 = 10, 5 = 30 


9 


2.0 x 10~ 3 


2.22 x 10~ 3 


1.5 x 10~ 3 


18 


2.26 x 10" 4 


2.66 x 10~ 4 


1.60 x 10~ 4 


36 


2.16 x 10" 5 


3.46 x 10~ 5 


1.32 x 10" 5 
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and we get in this case the results shown in table 5 and 6. 



Table 6: Maximum absolute errors for problem (|5.3p in yi. 



n 


\y(U) ~Vi\ 


10 


1.82 x 10" 1 


12 


2.33 x 10" 8 


15 


1.67 x 10" 8 



Conclusion 

In this paper we showed that hierarchical cascade models, motivated by turbulent flows, 
can lead to initial value problems of iV-th order, where N is the number of cascade steps. As 
an example we considered N = 7 and designed the optimum strategy to numerically treat 
the corresponding initial value problem. A non-polynomial spline method was developed for 
this. The numerical algorithm depends on some parameters a, /?, 7, 5 for which we derived 
explicit formulas. The method is observed to be second-order convergent for arbitrary 
choices of the parameters a, /?, 7 and 5 such that a + /? + 7 + <5 = 60 but if a, /?, 7 
and 5 are chosen as a = ^ — |, j3 = — + 5, 7 = — ^, then the method is of 
order five, due to the use of improved order end conditions. This is the optimum choice to 
numerically deal with cascade-like models of this type. In turbulent flows, driving forces 
and velocity differences often behave chaotically, corresponding to exponential growth for 
a limited amount of time. Consequently, we tested our method for these and other types 
of driving forces. We found that the maximum numerical errors observed are indeed very 
small if improved order end conditions are used. 
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